%% Author: Administrator
%% Created: 2009-2-27
%% Description: TODO: Add description to geo_grid
-module(geo_util).

-include("map.hrl").
%%
%% Include files
%%
-include_lib("eunit/include/eunit.hrl").

%%
%% Exported Functions
%%
-export([longitude_diff/1]).
-export([get_coordinates_from_address/1, get_quadtree_address/2, get_quadtree_address/3, cell_to_json/1]).

%%
%% API Functions
%%
%to_quardkey({Longitude, Latitude}) ->

longitude_diff(Latitude) ->
    Lat = Latitude * math:pi() / 180,
    111.413 * math:cos(Lat) - 0.094 * math:cos(3 * Lat).

sum_cell(Level) ->
    erlang:trunc((math:pow(4, Level) - 1) / 3).
%%
%% Google Maps tile code
%% see: http://intepid.com/stuff/gmkh/
%%
get_quadtree_address(Long, Lat) ->
    get_quadtree_address(Long, Lat, 24).

get_coordinates_from_address([$t|Id] = _CellId) ->
    coordinates_from_address(1, 0, 0, Id).

%%
%% get google maps image by: http://kh.google.com/kh?v=3&t={get_quadtree_address(Long, Lat, ZoomLevel)}
%% for example: http://kh.google.com/kh?v=3&t=trstrqqrtrqrrtqqq
%%
get_quadtree_address(Long, Lat, ZoomLevel) when ZoomLevel > 0 ->
    % now convert to normalized square coordinates
    % use standard equations to map into mercator projection
    X = (180.0 + Long) / 360.0,
    Y = mercator_to_normal(Lat),
    cat_address(ZoomLevel-1, X, Y, [$t]).
    
cell_to_json(Cell) ->
    Id = Cell#cell.id,
    CellZone = Cell#cell.zone,
    {struct, [
          {id, iolist_to_binary(Id)}, 
          {scale, Cell#cell.scale_factor},
          {zone, [CellZone#zone.left, CellZone#zone.bottom, CellZone#zone.right, CellZone#zone.top]}
    ]}.

%%
%% Local Functions
%%
coordinates_from_address(Scale, X, Y, []) ->
    MinLongitude = (X - 0.5) * 360.0,
    MinLatitude = normal_to_mercator(Y),
    MaxLongitude = (X + Scale - 0.5) * 360.0,
    MaxLatitude = normal_to_mercator(Y + Scale),
    _Longitude = (X + Scale * 0.5 - 0.5) * 360,
    _Latitude = normal_to_mercator(Y + Scale * 0.5),
    #zone{left=MinLongitude, right=MaxLongitude, top=MaxLatitude, bottom=MinLatitude};
coordinates_from_address(Scale, X, Y, [Char|RestId]) ->
    NewScale = Scale * 0.5,
    if 
        (Char == $r) or (Char == $s) ->
            NewX = X + NewScale;
        true ->
            NewX = X
    end,
    if 
        (Char == $t) or (Char == $s) ->
            NewY = Y + NewScale;
        true ->
            NewY = Y
    end,
    coordinates_from_address(NewScale, NewX, NewY, RestId).
    
normal_to_mercator(Y) ->
    Y3 = math:exp((Y - 0.5) * 2 * math:pi() * 2),
    Y4 = (Y3 - 1) / (Y3 + 1),
    math:asin(Y4)* (-180) / math:pi().
    
mercator_to_normal(Y) ->
    Y1 = math:sin(- math:pi() * Y / 180),
    Y2 = (1 + Y1) / (1 - Y1),
    0.5 + 0.5 * math:log(Y2) / (2 * math:pi()). % scale factor from radians to normalized 

cat_address(0, _X, _Y, Quard) ->
    lists:reverse(Quard);
cat_address(Digits, X, Y, Quard) ->
    X1 = X - erlang:trunc(X),
    Y1 = Y - erlang:trunc(Y),
    Ch = qrts(x_idx(X1) + y_idx(Y1)),
    cat_address(Digits-1, X1 * 2, Y1 * 2, [Ch|Quard]).

x_idx(X) when X >= 0.5 ->
    1;
x_idx(_X) ->
    0.

y_idx(Y) when Y >= 0.5 ->
    2;
y_idx(_Y) ->
    0.
    
qrts(0) ->
    $q;
qrts(1) ->
    $r;
qrts(2) ->
    $t;
qrts(3) ->
    $s.

normal_test() ->
    ?assertEqual(-90.0, normal_to_mercator(8)),
    ?assertEqual(-90.0, normal_to_mercator(19)),
    ?assertEqual(-89.98266599833198, normal_to_mercator(1.9)),
    ?assertEqual(85.05112877980656, normal_to_mercator(0)),
    ?assertEqual(89.98266599831095, normal_to_mercator(-0.9)).

mercator_test() ->
    ?assertEqual(0.3785791577410809, mercator_to_normal(40)),
    ?assertEqual(0.41257521185848506, mercator_to_normal(30)).

address_test() ->
    ?assertEqual("trstrqqrttsrsrstrqtsttss", get_quadtree_address(116, 39)),
    ?assertEqual("trstrqqrtrqrrtqqqqqqqqqq", get_quadtree_address(116.279296875, 39.87601941622491)).
